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Abstract —Recently, sum-of-squares (SOS) based methods 
have been used for the stability analysis and control synthesis 
of polynomial dynamical systems. This analysis framework was 
also extended to non-polynomial dynamical systems, including 
power systems, using an algebraic reformulation technique 
that recasts the system’s dynamics into a set of polynomial 
differential algebraic equations. Nevertheless, for large scale 
dynamical systems this method becomes inapplicable due to its 
computational complexity. For this reason we develop a subsys¬ 
tem based stability analysis approach using vector Lyapunov 
functions and introduce a parallel and scalable algorithm to 
infer the stability of the interconnected system with the help 
of the subsystem Lyapunov functions. Furthermore, we design 
adaptive and distributed control laws that guarantee asymptotic 
stability under a given external disturbance. Finally, we apply 
this algorithm for the stability analysis and control synthesis of 
a network preserving power system. 

I. INTRODUCTION 

Recently, a methodology for the algorithmic construction 
of Lyapunov functions for the transient stability analysis of 
classical power system models was introduced [1]. The pro¬ 
posed methodology uses advances in the theory of positive 
polynomials, semidefinite programming, and sum of squares 
decomposition, that have provided powerful nonlinear tools 
for the analysis of systems with polynomial vector fields 
[2]—[7]. In order to apply these techniques to power systems 
described by trigonometric nonlinearities an algebraic refor¬ 
mulation technique to recast the system’s dynamics into a 
set of polynomial differential algebraic equations was used 
[1], [4]. However, the sum-of-squares (SOS) approach does 
not scale well and, consequently, SOS methods only work 
for small systems with only a few state variables. For a large 
scale power system, it is then important to devise a subsystem 
based stability analysis approach using the concept of vec¬ 
tor Lyapunov functions and the decomposition-aggregation 
method [8], [9]. 

Formulations using vector Lyapunov functions [10], [11] 
are computationally attractive because of their parallel struc¬ 
ture and scalability. In [9], it was shown that if the subsystem 
Lyapunov functions and the interactions satisfy certain con¬ 
ditions, then application of comparison equations can provide 
a certificate of exponential stability of the interconnected 
systems. The approach we propose here uses instead a 
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generalization to interconnected systems of SOS based robust 
stability analysis techniques developed to estimate the effect 
of parametric uncertainty [5] and external disturbances on a 
polynomial system [7]. 

In this work we use sum-of-squares analysis methods to 
devise a new algorithmic certification of asymptotic stabil¬ 
ity via the vector Lyapunov function approach. While this 
approach is generic, we apply the proposed algorithm to 
analyze the stability of a structure preserving power system 
model [12], [13]. Unlike the network preserving models stud¬ 
ied in the literature we do allow for nonzero transfer conduc¬ 
tances in the transmission lines. The network is decomposed 
into a number of low order interacting subsystems. For each 
such subsystem, a SOS based expanding interior algorithm 
[1], [7] is used to obtain estimate of region of attraction as 
sub-level sets of polynomial Lyapunov functions. Finally a 
sum-of-squares based scalable and parallel algorithm is used 
to certify stability in the sense of Lyapunov of the intercon¬ 
nected system by using the subsystem Lyapunov functions 
computed in the previous step. A distributed control strategy 
is proposed that can guarantee asymptotic stability of the 
interconnected system under given disturbances. 

Following some brief background in Sec. II we outline the 
problem statement in Sec. III. An algorithmic approach to 
certifying asymptotic stability is presented in Sec. IV while 
a distributed control strategy is discussed in Sec. V. Sec. VI 
shows an application of our stability analysis and control 
approach to a network of three generators and six frequency 
dependent loads. We conclude the article in Sec. VII. 

II. BASIC CONCEPTS AND BACKGROUND 

We will first discuss how the stability of a dynamical 
system can be analyzed by constructing suitable Lyapunov 
functions. Then we briefly refer to sum-of-square polynomi¬ 
als and a very useful result which helps us in formulating 
the sum-of-squares problems. 

A. Lyapunov Stability Methods 

Let us consider the dynamical system described by the 
following polynomial differential algebraic equations (DAE) 

z = F(z), (la) 

0 = G(z), (lb) 

where z G M m , and F : M m -A M m , G : M m -A M. q are vectors of 
polynomial functions. We assume without loss of generality 
that the origin is a stable equilibrium point for this system, 
i.e. F(0) =0 and G(0) = 0. 



The following extension of Lyapunov stability theorem 
to dynamical systems described by differential algebraic 
equations presents a sufficient condition of stability through 
the construction of a certain positive definite function [1], 

[5]. 

Theorem 1 If there exists an open set D C M m , with D a 
semi-algehraic domain defined by the following inequality 
and equality constraints, 

D = {z G M m | P - p{z) > 0,G(z) = 0}, (2) 

with p(z) a positive definite polynomial and > 0 to ensure 
that D is connected and contains z = 0, and a continuously 
differentiable function V : D —M such that V (0) = 0, and 

V(z)>0,VzGD\{ 0}, (3a) 

-V(z)>0,VzGD\{ 0}, (3b) 

where V(z) = VV r • F(z), then z = 0 is an asymptotically 
stable equilibrium of ( 1 ). 

When there exists such a function V(z) 9 the region of 
attraction (ROA) of the stable equilibrium point at origin 
can be (conservatively) estimated as 

^:={zG/)|y(z)<r} (4a) 

where, f f nax argmax {z G D \V (z) < 7 } C D (4b) 
7 

Without any loss of generality, the Lyapunov function can 
be scaled by y max 9 so that the ROA is given by, 

^ A :={zeDlV(z)<lj (5) 

Henceforth, we would assume that the ROA is estimated to 
be the sub-level set of V (z) = L 

B. Sum-of-Squares and Positivestellensatz Theorem 

While Thm. 1 gives a sufficient condition for stability, it is 
not a trivial task to find a function V (z) that satisfies the con¬ 
ditions of stability. Relatively recent studies have explored 
how sum-of-squares (SOS) based optimization techniques 
can be utilized in finding Lyapunov functions by restricting 
the search space to SOS polynomials [1], [2], [7], [14]. Let 
us denote as the set of all polynomials in z G M m . Then, 

Definition 1 A multivariate polynomial p(z) G S% m is a sum- 
of-squares (SOS) if there exist some polynomial functions 
Hi (z), i == 1 .. .r such that p(z) = JT = 1 //^(z), and the set of 
all such SOS polynomials is denoted by 

£m ■= {P(Z ) G | P IS SOS } . ( 6 ) 

Given a polynomial p G & m , checking if it is SOS is a semi- 
definite problem which can be solved with a MATLAB® 
toolbox SOSTOOLS [15], [16] along with a semidefinite 
programming solver such as SeDuMi [17]. 

Hence, if we search for a polynomial Lyapunov function in 
Thm. 1, and if we relax the polynomial non-negativity con¬ 
ditions to appropriate polynomial sum of squares (SOS) con¬ 
ditions, testing SOS conditions can then be done efficiently 


using semidefinite programming (SDP) [2]. Moreover, an 
important result from algebraic geometry, called Putinar’s 
Positivstellensatz theorem [18], [19], helps in translating the 
SOS conditions into SOS feasibility problems. Before stating 
the theorem, let us define: 

Definition 2 Given gj G M m , for j = l,2,...,r, the 
quadratic module generated by gj’s is ^(gi,g2,• • • ,gr) : = 

{oo+Ey-=i Ojgj |oo, oj e r m ,V;} • 

Then the Putinar’s Positivestellensatz theorem states [19], 

Theorem 2 Let J(f — {z G W 1 |gi(z) > 0,.. .,g r (z) > 0} be 
a compact set. Suppose there exists u(z) G such that 

u{z) £^(ghg 2 ,---,gr), (7a) 

and, {z G W 1 \ u(z) >0} is compact. (7b) 

If p(z) is positive on J(f, then p(z) G ^(gi,g2? • • • ,gr)- 

It can be shown that for all gf s used in this work the 
constraints (7) would be redundant, i.e. the existence of u(z) 
would be guaranteed. If we further note that the equality con¬ 
straints defined by the components of G(z) can be expressed 
by the pairs of inequalities gk > 0 , gk < 0 , k = 1 ,..., q, and if 
we define g q+ \ = fi — p, then using Thm. 2, with r = q-\- 1, 
the search for V (z) becomes a search for a feasible solution 
of the following problem with SOS constraints [1], [7]: 

V-^G-s 1 (fi-p)-<ln€l m ( 8 a) 

-V-XlG-S2^-p)-<t> 2 ei: m ( 8 b) 

where st G Z m , G are vectors of polynomial functions, 
and we choose tyfiz) = £iT!k=i z b £ i > * = 

Similarly, the search for y 10 * in Eq. (4) can be formulated 
as a SOS programming problem and solved using a bisection 
search over 7 — see [1], [7] for details. 

III. PROBLEM OUTLINE 

Given the full dynamical system (1), we seek to decom¬ 
pose it into S weakly interacting subsystems as 

Zi = F i (z i )+H i {z), z=l,2,...,S, (9a) 

0 = G t (zi ), (9b) 

Fi(0)m Hi( 0)=0 (9c) 

where Zi G M m/ is the state of the i -th subsystem Si, Ff s 
denote the isolated subsystem dynamics, G ; - G are vec¬ 
tors of polynomials defining the algebraic constraints of the 
subsystem, and Hf s are the interactions from the neighbors 
[20]. If we denote by 3?i the subset of state points that 
correspond to z,i then their union, i.e. 

= jFiUjr 2 u...u jg, (io) 

where 3? = {zuZ2, • • • represents the set of state points 
of the whole dynamical system (1). In the decomposition 
used in this paper the subsets 3%, i = 1,..., S, are not disjoint. 





We assume that the interactions can be expressed as, 

v/. //,(.-) m£ H ij(zi,Zj) (11) 

j^i 

where H[j quantifies how the states zj of subsystem j affect 
the dynamics of Zi . Let us also denote by 

9% ■= (0 U {y | 3 Zi,Zj, s.t. Hij ( Zi,Zj ) ^ 0 } , (12) 

the set of neighbors of node i (including the subsystem i 
itself). We assume that the isolated subsystems are individu¬ 
ally (locally) stable, and there exist Lyapunov functions for 
each of the isolated subsystems. The goal is to develop a 
framework for the stability analysis of the full interconnected 
system by using the local subsystem Lyapunov functions 
and considering the neighbor interactions. In the work re¬ 
ported here we use an ’’extreme” decomposition in which 
the subsystems are defined by the nodes of the original 
network together with a reference node that is shared by 
all subsystems. A similar approach has been used in [21]. 
A decomposition algorithm proposed in [22], and used for a 
classical power system model in [23], was not used here, 
but will be used in the future to analyze the impact the 
decomposition has on the performance of the algorithm. 

Given a decomposition (9), our next goal is to find poly¬ 
nomial Lyapunov functions Vfn) G M mi for each isolated 
subsystem i = 1,2,....5, 

Zi = Fi(zi ), (13a) 

0 = Gi( Zi ). (13b) 

The first step in this search is to solve the SOS program ( 8 ) 
which for each subsystem -9 is formulated as 

Vi - Gi - s n (A ~Pi)~ Ai e E m; (14a) 

-Vi - A/(G, - sa(fii ~Pi) ~ te e E mi (14b) 

where A > 0, sn,sa G £ m; , An,A ;2 G and 

0,1 (zi), ( pa(zi ), Pi{zi ) are positive definite polynomials. 
Starting from an initial Lyapunov function candidate 
obtained by solving (14), and a corresponding estimate 
of the region of attraction, an iterative process called 
expanding interior algorithm , [1], [7], is used to iteratively 
enlarge the estimate of the region of attraction by finding 
a better Lyapunov function at each step of the algorithm. 
At the completion of this iterative step, the stability 
of each isolated subsystem (assuming no interaction) 
is quantified by its Lyapunov function V/(z/), with an 
estimation of the boundary of the domain of attraction given 
by ^A,« = {zi G W* | Gi(zi) = 0 ,Vi(zi) < 1}. 

The Lyapunov level-sets can be used to express the 
strength of a disturbance. The equilibrium point of the system 
at origin corresponds to the level set V/(0) = 0, V/. If there is 
a disturbance from this equilibrium point, the states of the 
system would move to some point x( 0 ) away from the origin. 
This disturbed initial condition would result in positive level- 
sets Vi(zi( 0)) = tP G (0,1] for some or all of the subsystems. 


A necessary and sufficient condition of asymptotic stability 
can then be translated into the condition 

Vi, v;-(zi(0)) = => Vi, lim Vi(zi(t)) = 0 (15) 

t —»+oo 

In the rest of the article, we present SOS algorithms to 
test stability conditions and design local (subsystem-level) 
control laws to achieve asymptotic stability. 

IV. STABILITY UNDER INTERACTIONS 

It is assumed that the isolated subsystems in (13) are all 
(locally) asymptotically stable, and that there exist subsystem 
Lyapunov functions V/(z/), V/. The estimated region of attrac¬ 
tion of the interconnected system under no interaction, is 
given by the cross-product of the regions of attraction of the 
isolated subsystems, which are defined as sub-unity- 
level sets of the corresponding (properly scaled) subsystem 
Lyapunov functions (as in (4)-(5)), i.e. 

:= &AA X &A2 X • • • X ,m 

where, &A,i=Ui G M m/ |G,-(z/) = 0,V/(z,-) < 1}, VL 

In presence of non-zero interactions, the resulting ROA 
would be different. If there exists a Lyapunov function for 
the interconnected system, the ROA for the whole system 
could be expressed as some sub-level set of that Lyapunov 
function. While it is very hard to obtain a scalar Lyapunov 
function for the full interconnected system, one could use 
vector Lyapunov function approach to obtain certification of 
stability in a scalable way. 

In this present work, we choose not to impose any further 
restriction on the Lyapunov functions Vfzf) than requiring 
that those are in polynomial forms, and concern ourselves 
with asymptotic stability. We now present a distributed 
iterative procedure which can be used to certify asymptotic 
stability in a domain defined by sub-level sets of the subsys¬ 
tem Lyapunov functions. 

A. Algorithmic Test of Asymptotic Stabiltiy 

Before proceeding to explaining our algorithm, let us first 
note the following result: 

Lemma 1 Suppose, for all i G {1,2,... ,5}, there exists 
a strictly monotonically decreasing sequence of scalars 
{ef } , k G {0,1,2,...}, such that 

Vi,*, Vi(zi):=VVi(zi) T (Fi(zi)+Hi(z))<0 , VzG^f (17a) 

{ £ k + l < Vi(zi) < £^, } 

z G IT'' Vj(zj) < £j Vy G A\ {*} > (17b) 
Gj{ Z j) = 0 Vy G jV[ J 

Then the system (1) is asymptotically stable in the domain 
{z G M m | nf=i Vi(zi) < ef }, if l im^ +0 o£f = 0,V/ If the limit 
condition does not hold, we can only guarantee stability in 
the sense of Lyapunov [24], [25]. 

Proof: Please refer to Appendix A. ■ 




Using the Lemma 1 we can devise a simple iterative SOS 
algorithm to certify whether or not a domain $ defined by 


0 := 



m 

n{^fc)<^,G J -fe)=o} 


( 18 ) 


system. Moreover, the algorithm motivates the design of 
a distributed control strategy that can ascertain asymptotic 
stability. 


V. LOCAL CONTROL 


for some scalars £ (0,1], V/, is a region of asymptotic 
stability for the system in (1). It is to be noted that, using 
the Putinar’s Positivstellensatz theorem (Theorem 2), the 
condition in (17) essentially translates into equivalent SOS 
feasibility conditions 

Vi,*, s.t. 

- vvf OR + Hi) -of (v;--ef +1 ) 

- E ( e i -- E e < 19 ) 

jeA jeA 

where, m ; -. 

The algorithmic steps to ascertain asymptotic stability are as 
outlined below: 

1) We initialize £? = yf,V/ £ {1,2,... ,5}, and choose a 
sufficiently small £ £ M + . 

2) At the start of the k-th iteration loop, we assume to 
know the scalars {£?,£/,... ,£^} ,Vz, and our aim is 
to compute the scalars £* +1 , Vi such that (19) holds. 
Essentially we want to solve the optimization problem, 

Vi, min £f +1 s.t. (19) holds. (20) 


In this section, we discuss design of a local and minimal 
control strategy such that the system in ( 1 ) is asymptotically 
stable in a domain Q) defined in (18). We use the term 
minimal to suggest that the control be applied only in certain 
regions, and not everywhere, in the state space, while by the 
term local we suggest that the control is computable and 
implementable on a subsystem level. 

We envision the control to be computed by each subsystem 
at each iteration loop. At k-th iteration, \/k £ {0,1,2,...}, 
the i-th subsystem, Vi £ {1,2,... ,5} performs the following 
tasks: 

1) It identifies if it belongs to the following set 



/ 

W? {Fi + Gi) > 0, ) 

W k := < 


Vi = & \ 

Vj < e) Vj G jVi\ {/} 



G ; '(.-/)=.()V/C. (,• J 


which can be checked locally. If i ^ control is not 
necessary and it sets Uf = 0 and proceeds to task 3. 
If, however, i £ % k , it proceeds to task 2 to compute 
a control law. 

2) If i £ % k , a polynomial state-feedback control law U\ : 
M m/ -£ M m % U\( 0) = 0, is computed such that 


This is solved by performing a bisection search for 
minimum £^ +1 over the range [ 0 ,£*]. 

If (20) is infeasible at 0-th iteration for any i £ 
{1,2,... ,5}, we conclude that the system cannot be 
guaranteed to be asymptotically stable in Q), and abort 
the iteration. Otherwise we move on to step 3. 

3) If (z\ — £f +1 ) > £, V/, we continue from step 2 for the 
(k+l)-th iteration loop. Otherwise we stop the iteration 
deciding that the limits of the sequences \z\ } , Vi, have 
been attained. Further, if the limits are all zero, we 
certify asymptotic stability in Q). 

B. Remarks 

The algorithm presented in Sec. IV-A describes how one 
can determine asymptotic stability of an interconnected sys¬ 
tem in a domain Q) defined by the subsystem sub-level sets. 
This test can be performed locally, and in a parallel way, at 
each subsystem level. The Lyapunov functions V/’s are to be 
found before the start of the analysis, and communicated to 
the neighboring subsystems. Then during each analysis, it is 
assumed that the neighboring subsystems can communicate 
with each other the computed sequences \z \} in real-time. 
With the help of the stored Lyapunov functions, and the 
updated {e k } of the neighbors, each subsystem will continue 
the iterative process outlined in Sec. IV-A. Since only the 
neighbor information is required, this algorithm is reasonably 
scalable with respect to the size of the full interconnected 


vvf (Fi+Gi + U k ) < 0, Vz € +, 


where f k 



Vi = 

Vj<e+je+\{i} 
Gj( Zj ) = 0Vje^ 


( 21 ) 


This produces the equivalent SOS condition, 

- VVf (Fi + Gi + uf) - p\ (ef - v;-) (22) 

- E °y (< £ j ~ V J ) - E +Y G i e 

jeA\{i] jeA 

Pi G i j £ V/' 7^ i■ hj G , ifli = ^ fflj 

jeA 


3) Finally, it performs the search over minimum £* +1 , 
as in ( 20 ), with the un-controlled subsystem dynamics 
(Fi + Gi) in the feasibility condition (19) replaced by 
the controlled dynamics (F t + Gi + U k ). 

To summarize, each subsystem i computes control laws 
jjk . ^Mi jjk j'q) _ (tunng each £_th iteration, 

so that the subsystem dynamics under control becomes: 


Vi, V*, Vx G Sf, 


J FiM+Hifa), Gi(zi)= 0, 

\ Fi(zi)+Hi(zi) + Ut(zi), Gi{zi)= 0, 


i <£ <fy k 
i G ^ k 


where were defined in (17b). 



A. Remarks 

Often it is important to impose certain additional con¬ 
straints on the possible control laws, such as bounds on 
the control effort. Although control bounds can be easily 
incorporated in the SOS formulation, we decide to keep that 
for future studies. However, that since we apply controls Uf 
only on certain subsystems i G and in certain domains 
@)\ C 3), the control effort would be reasonably bounded. 

VI. RESULTS 

Let us describe the model of the interconnected system that 
we use here, and two examples to illustrate the applications 
of the stability analysis algorithm and control design. 

A. Model Description 

We will consider the network preserving model with linear 
frequency dependent real power loads introduced in [ 12 ], 
[13]. The network consists of G generators and L load buses 
connected by transmission lines. We number the load buses 
1 , 2 ,... ,L and the generator buses L+ 1 ,...where n = 
L + G is the total number of nodes. The node voltages are 
denoted by E\ Z8 \,..., E n Z8 n , where <5i,..., 8 n are the phase 
angles and the magnitudes E\,...,E n are assumed constant. 
The electrical power injected into network at node i is 

P Ei (S) = EfGu + £ EiEjYij cos(5,- - Sj - %) (23) 

jJAi 

where Uy is the modulus, and 6jj the phase angle, of the 
transfer admittance between nodes i and j. 

Lor small frequency variations around the operating point 
Pdi the dynamics of the load nodes are described by 

DiSi = -P Di -P Ei {S ), i = 1,... ,L, (24) 

where Di > 0 is the load-frequency coefficient and Po i > 0 
is the real power drawn at the load buses. Each generator 
dynamics are modeled by the swing equations. Thus, for 
i = 1 ,... ,G, 

^L+f = tO/, (25 a) 

Mj (bj + D E +i (Di = Pm l — PE L+i ( 8 ), (25b) 

where > 0 is the generator inertia constant, D E +i > 0 is 
the generator damping coefficient, Pm, > 0 is the mechanical 
power input. 

We assume that the dynamical system has a stable equi¬ 
librium point given by (S s ,(O s = 0 ) where 8 S is the solution 
of the following set of nonlinear equations, 

-P Dl ~P Ei {8 s ) = 0,for i = 1,... ,L (26) 

Pm, ~P EL+i ($s ) = 0,for i = 1,... ,G (27) 

The state space of the dynamical system is described by 
the relative angles 8f n = 8j■ — 8 n , for i = 1,..., L + G — 1, with 
respect to a reference node n (generator G, considered to 
have the largest inertia), and relative generator speeds (Di n = 
CQj — co n , for i = 1,..., G — 1, for all generators except the 
reference generator for which we consider the absolute speed 
(O n = (Dg- (We have considered that the ratio Di+i/Mi is 


uniform.) The dynamics of the relative angles and speeds 
are obvious and are not explicitly presented here. Linally, 
we make the following change of variables, 8 8 + 5?, in 

(24) and (25) in order to transfer the stable equilibrium point 
to the origin in phase space. 

B. Recasting the Power System Dynamics 

SOS programming methods cannot be directly applied to 
study the stability of power system models because their 
dynamics contain trigonometric nonlinearities and are not 
polynomial. Lor this reason a systematic methodology to 
recast their dynamics into a polynomial form is necessary [3], 
[4]. The recasting introduces a number of equality constraints 
restricting the states to a manifold having the original state 
dimension. Lor the network preserving power system model 
recasting is achieved by a non-linear change of variables, 

<*■>->{ S= i-tS)} <2M 

{(Din} -A {z2(L+G—\)+i = ^in} , i = 1, •. •, G — 1 , (28b) 

{&>;} —>> {Z2(L+G-1)+; = G = G, (28c) 

and the introduction of the following constraints 

4-1 + 4, - 222, = o, / = 1,..., L + G - 1. (29) 

Thus, recasting produces a dynamical system with a larger 
state dimension, z G M m , where m = 2(L + G— 1) + G and 
introduces q = L + G — 1 equality constraints. The stable 
equilibrium point of the original system is mapped to = 0. 
The original system dynamics are recasted into the polyno¬ 
mial differential algebraic equations (1). 

C. Test Case: IEEE 9-bus System 



Fig. 1. Network of three generators and six load nodes. We perform an 
overlapping decomposition in which the speed dynamics of the reference 
node (generator node 3) is shared with all the subsystems. This results in the 
following overlapping subsystems: S\ = (4,1 },^2 = (5,1}A3 = (6,1},^4 = 
(7, 1},S 5 = (8, 1} ,S 6 = {9, 1} A 7 = {2, 1},S 8 = (3, 1}. 

We will be using the Western System Coordinating Coun¬ 
cil (WSCC) 9-bus system (commonly, the IEEE 9-bus model 
[26]), for our analysis. To better illustrate the scope of 
our approach, we redistribute the loads so that each of the 
nodes 4, 6 and 9 in the network has certain loads, as shown 
















(a) ROA of subsystem 2 (node 5) (b) ROA of subsystem 7 (node 2) 

Fig. 2. Estimates of regions of attraction of isolated subsystems using expanding interior algorithm. 





(a) Load angles, before control 


(b) Generator states, before control 


(c) Subsystem Lyapunov functions 


Fig. 3. System states and Lyapunov functions under a disturbance, before control. 


TABLE I 

Algorithm certifies asymptotic stability, under control 


k 

¥ 

7 £ 

£ 2 

£ 3 

p* 

fc 4 

pfc 

b 6 

¥ 

71 

£ 8 

0 

0.0001* 

0.0007* 

0.0002* 

0.3471 

0.0003* 

0.0001* 

0.6663* 

0.0001* 

1 

0.0000 

0.0000 

0.0002 

0.0312 

0.0000 

0.0001 

0.4451 

0.0000 

2 

0.0000 

0.0000 

0.0001 

0.0203 

0.0000 

0.0000 

0.4260 

0.0001 

3 

0.0000 

0.0000 

0.0001 

0.0195 

0.0000 

0.0000 

0.3478 

0.0000 

4 

0.0000 

0.0000 

0.0000 

0.0154 

0.0000 

0.0000 

0.3383 

0.0000 










11 

0.0000 

0.0000 

0.0000 

0.0103 

0.0000 

0.0000 

0.2481 

0.0000 

12 

0.0000 

0.0000 

0.0000 

0.0010 

0.0000 

0.0000 

0.2374 

0.0000 

13 

0.0000 

0.0000 

0.0000 

0.0001 

0.0000 

0.0000 

0.0492 

0.0000 

14 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0094 

0.0000 

15 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 



t[s] 



t[s] 



t [S] 


(a) Load angles, with control at Gen. 2 (b) Generator states, with control at Gen. 2 


(c) Subsystem Lyapunov functions 


Fig. 4. System states and Lyapunov functions under local and temporary control applied at generator 2, for the duration t G [3s, 3.127s]. 





























































































in Fig. 1. Thus, in our modified network, each node has 
some associated dynamics. Further, using an overlapping 
decomposition [27], we construct eight subsystems where 
each subsystem is composed of either a load node or a (non¬ 
reference) generator node, along with the reference generator 
node. Essentially this requires that the reference generator 
speed value is to be communicated to all the subsystems, 
in real-time 1 . We choose D^yi = 5M ? for the generators and 
select Dj randomly from [1,2] for the loads. 

We use the expanding interior algorithm to find estimates 
of the ROA for each isolated subsystem. Fig. 2 shows the 
evolution of ROA estimates projected on co n = 0. In Fig. 2(a) 
we see how the level set of the Lyapunov function for a load 
node evolves as the ROA estimate, computed as in (4)-(5), 
grows. Fig. 2(b) shows a comparison of the true ROA for a 
generator node with a sequence of its estimates where the 
outermost ‘blue’ contour represents the final estimate 2 . 

Any randomly picked initial condition, jt(0) G where 
is defined in (16), can be mapped into corresponding 
subsystem Lyapunov function level sets, yP = Vi(zi(0)),Vi. 
Then, by choosing ef = yf , we apply the iterative stability 
analysis algorithm to determine whether or not the domain Q) 
in (18) is a region of asymptotic stability, and if not, compute 
the necessary control by (21). 

D. Disturbance Analysis 

We assume that initially the network was operating at 
equilibrium. A disturbance was created by tripping the line 
between nodes 5-7 for the duration t G [0, 3s] and also 
tripping the line between nodes 7-8 for t G [Is, 3s], essentially 
disconnecting the nodes 7 and 2 from the rest of the network 
for t G [Is, 3s]. The end of the disturbance provides the initial 
condition for the stability study. The evolution of states and 
the Lyapunov functions from this initial condition (i.e., at t = 
3s) is shown in Fig. 3. The initial condition is asymptotically 
stable. Further, the local effect of the disturbance is clear 
from Fig. 3(c), where apart from V^(t) and Vj(t) 9 all other 
Lyapunov functions are very small. 

In Table I the results from an application of the stability 
analysis algorithm to the initial conditions are shown. The 
row corresponding to k = 0 lists the initial level sets, while 
subsequent rows list the sequence {£*} for all the subsys¬ 
tems. The * shows when the algorithm feels the necessity of 
applying control 3 to guarantee the decrease of level set. For 
example, at the first iteration step, the algorithm prescribes 
applying state-feedback control 

^' 7 ° = 5.1cos52,i-38.4ft>2 i i-1.8ffli-71.4sin52,i-5.1 (30) 

which is applied to the speed dynamics of the generator 2 
of Sj. It can be seen from Fig. 3(c) that initially Vj(t) shows 
a slight increase before starting to decrease, which is why 

Communication bandwidth and time-delays associated with such a 
model are important issues that are beyond the scope of the present work. 

2 The ROA estimates are obtained using quadratic polynomials. 

3 We chose to seek only linear controllers in the z variables which are 
applied on the angle dynamics for the loads, and on the speed dynamics for 
the generators. The controllers are nonlinear in the original state space. 


the control is prescribed. For the same reason, control is also 
suggested for subsystems other than 54 . However, other than 
5*7, control action can be safely ignored because the level 
sets are very close to zero. 

After k = 1, no control is deemed necessary by the 
algorithm, and finally after iteration k = 15, the certificate of 
asymptotic stability is obtained. Fig. 4 shows the evolution 
of states and Lyapunov functions under the action of control 
applied at the generator 2. The control is applied only for the 
time during which 0.4451 < Vj(t) < 0.6663, which amounts 
to the time interval t G [3s, 3.127s], as shown in Fig. 4(c). It 
is to be noted that under the control action, Vj(t) becomes 
negative at t = 3s. 

VII. CONCLUSIONS 

In this work, we present a distributed algorithmic approach 
to infer the stability of an interconnected power system 
under a given disturbance. When this stability cannot be 
certified we design local and minimal control laws that 
guarantee asymptotic stability. The approach presented here 
is parallel and scalable. While this method is applied here 
considering dynamic loads in a network preserving power 
system model, it can be extended to more general, possibly 
more complex, dynamical representation of the power system 
network. Future work needs to address the issues of bounded 
control effort, voltage dynamics of loads and inclusion of 
available control mechanisms, such as speed governors. 

Appendix 


A. Proof of Lemma 1 

We note that since lim^ +00 £f = 0,V/, 


V<5 G 0,min£- 


, 3K, s.t. ef < 8 Vfc > Kyi. (31) 


Let us assume, without any loss of generality, that 

s 


n{vK*0<s p > <?.•(*) so} 

i=l 


3to > 0, S.t. z(*o) G jzG 

Then, 

Vi, Vi(t) < ef + f Vi(r)dT, Vt > t 0 
Jt 0 

=> 3 tf < t () + (ef - ef) / rj, rj := sup V (x) < 0 

x<E@f 

s.t. Vi(t) < ef, Vt >tf . 

Hence, we can argue that, 

Vi{t) < ef, Vt > t 0 ,Vi 

==> 3 f 1 := maxf/, s.t. V,(t) < e}. Vf > tfVi. 

i 

Following similar arguments it is easy to show that, 

Vi(t) < ef, Vt > to, Vi 

=*> Vk, 3 t k > t 0 , s.t. Vi(t) < ef, Vt > t k , Vi. (32) 
Finally combining (31) and (32) we observe, 


V<5 G 0,min£? 


3 t K > to, s.t. Vi{t) <5, \/t > t K , VL 


which concludes the proof, because of (15). 
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